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Abstract 

The limits of reversible deformation in graphene under various loadings are examined using lattice-dynamical 
stability analysis. This information is then used to construct a comprehensive lattice-stability limit surface 
for graphene, which provides an analytical description of incipient lattice instabilities of all kinds, for ar¬ 
bitrary deformations, parametrized in terms of symmetry-invariants of strain/stress. Symmetry-invariants 
allow obtaining an accurate parametrization with a minimal number of coefficients. Based on this limit sur¬ 
face, we deduce a general continuum criterion for the onset of all kinds of lattice-stabilities in graphene: an 
instability appears when the magnitude of the deviatoric strain 7 reaches a critical value 7 *^ which depends 
upon the mean hydrostatic strain £ and the directionality 0 of the principal deviatoric stretch with respect 
to reference lattice orientation. We also distinguish between the distinct regions of the limit surface that 
correspond to fundamentally-different mechanisms of lattice instabilities in graphene, such as structural 
versus material instabilities, and long-wave (elastic) versus short-wave instabilities. Utility of this limit 
surface is demonstrated in assessment of incipient failures in defect-free graphene via its implementation 
in a continuum Finite Elements Analysis (FEA). The resulting scheme enables on-the-fiy assessments of 
not only the macroscopic conditions (e.g., load; deflection) but also the microscopic conditions (e.g., local 
stress/strain; spatial location, temporal proximity, and nature of incipient lattice instability) at which an 
instability occurs in a defect-free graphene sheet subjected to an arbitrary loading condition. 

Keywords: Graphene, ideal strength, lattice-stability limits, finite element analysis. 


1. Introduction 

Limits to reversible deformation—the stress or strain at which elastic-to-inelastic transition takes place 
in a material — pose fundamental constraints on a material’s performance and determine its strength. The 
limiting stress or strain, in general, depends upon the loading path, and the dependence is often described 
by means of a phenomenological model termed a limiting (or failure) criterion. Mathematically, a limiting 
criterion is represented as a surface in stress or strain space, which separates the stable states of reversible 
deformation from the Tailed’ or irreversibly-deformed states. For many materials, the limiting criterion is 
simple enough to be characterized by one or two material constants, which are readily determined from 
experiments. Examples include the Mises (or Tresca) yield criterion for metals, the Mohr-Coulomb criterion 
for cohesive-frictional solids m, the Drucker-Prager criterion for pressure-dependent solids and the 
Hoek-Brown criterion for rocks m- 

In crystalline materials that are free from lattice imperfections, the limit to elastic deformation sets an 
upper bound on the material strength [39]. This upper bound, termed the ideal strength, depends on the 
intrinsic nature of bonding between atoms in the material. Because of the various types of lattice defects, 
such as dislocations, grain boundaries, interstitial impurities, and voids which normally exist in materials 
uniiso], most conventional materials are irreversibly-deformed at stress-levels well below the ideal strength. 
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However, in recent years, advancements in nanotechnology have enabled fabrication and growth of defect- 
free two-dimensional crystals in which mechanical failure indeed occurs upon reaching stress levels near the 
ideal strength [HUM]. 

Graphene, an atomic monolayer comprising a hexagonal network of covalently-bonded C-atoms, is a 
representative example of such materials (see Fig. 0 ). Experimental studies have shown that defect-free, 
single-crystalline graphene can sustain near-ideal-strength stresses while remaining within the reversible 
regime of deforequationmation [321 Hi- Beyond the limit of elastic deformation, the fate of the material is 
determined by a strength-limiting mechanism such as incipient plasticity or crack-initiation. Since, in the 
absence of lattice-imperfections, a strength-limiting mechanism can only be activated by a lattice-instability 
[33l , the incipient failure of a defect-free crystalline material is intrinsically related to the loss of internal 

lattice-stability. The point at which the loss of lattice-stability occurs is called the lattice-stability limits and 
it varies from loading path to path. Further, the dependence of the lattice-stability limit on the loading 
condition in crystalline materials is generally too complex to be adequately characterized by one or two 
parameters. 

Individually, lattice-instabilities of all kinds can be assessed via the lattice-dynamical stability analysis 
(see Born & Huang [7]; Hill m), which asserts that a necessary and sufficient criterion for an ideal crystal 
under arbitrary uniform loading to be stable is that it exhibits stability with respect to bounded perturba¬ 
tions of all wavelengths. Integration of the lattice-dynamical stability analysis with a continuum analysis 
scheme such as FEA would be ideal for failure analysis of defect-free graphene crystals. This could enable 
assessment of incipient failures and ideal strength of graphene, under arbitrary loading conditions at realistic 
length-scales and slow enough loading rates, directly in a continuum-level simulation. However, stability 
analysis based on lattice-dynamics requires carrying out an elaborate sequence of computationally-expensive 
steps that can not be treated within the confines of an analytical framework, making it difficult to integrate 
the lattice-dynamical stability analysis into a continuum scheme. Therefore, there remains a need for a 
general continuum criterion that could describe the onset of instability (of any kind) in graphene under an 
arbitrary state of deformation. 

The aim of this work is to construct a comprehensive lattice-stability limit surface, which constitutes an 
analytical parametrization of incipient lattice instabilities of all kinds, over the space of all homogeneous 
deformations, in terms of stress/strain. There are two main difficulties in obtaining such a parametriza¬ 
tion: First, crystalline materials are intrinsically anisotropic, so material response, including lattice-stability 
limit, varies with orientation [33]. Secondly, two fundamentally-different types of lattice-instabilities govern 
strength-limiting mechanisms under different loading conditions [Ml uni: a long-wave (or elastic) instability 
and a short-wave (or soft mode) instability. The condition for onset of an elastic instability can be param¬ 
eterized in terms of strain via acoustic tensor analysis (see Kumar & Parks [25] for details), whereas the 
short-wave instabilities are much more complex since there is no continuum framework for parametrization 
of limiting conditions governing the onset of short-wave instabilities. 

The proposed parametrization, in order to overcome the above-mentioned difficulties, employs inter¬ 
polation of the lattice-stability limits of graphene, corresponding to some representative homogeneous de¬ 
formation modes, in the basis of symmetry-invariants of the strain/stress measure. Symmetry-invariants 
are special scalar functions of strain/stress that remain invariant under the point group operations. The 
use of symmetry-invariants ensures that the parametrization possesses the appropriate material anisotropy 
of graphene. It also introduces substantial functional simplification, reducing the requisite representative 
deformation modes needed for interpolation to a small set of biaxial deformations along the two special 
symmetry-directions in graphene: armchair and zigzag. The individual lattice-stability limits used in the 
interpolation are obtained by atomistic-level lattice-dynamical stability analysis based on phonons to ensure 
that lattice instabilities of all kinds are captured in the analysis. 

This paper is organized as follows. In Sec. 0. we briefly summarize the kinematics of graphene, outline 
the general framework of invariant-based representation theory of anisotropic materials, explicitly derive 
symmetry-invariants of the strain measure with respect to the symmetry-group of graphene, and discuss the 
implementation of these ideas in the context of an (incipient) instability function of graphene. Employing 
a representation of the ideas outlined in Sec. 0, we systematically determine an analytical form for the 
limit surface for graphene in Sec. 0 . Then, in Sec. 0 . we map the instability function from strain space 
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to stress space. The model is validated in Sec. in Sec.( |5.1|), w e assess the numerical accuracy of the 
representation of failure model based on interpolation; in Sec.( [5.2[ ), we compare the material instabilities 
predicted by the failure model with those from acoustic tensor analysis; and in Sec.(5.3), we compare the 
model predictions for buckling instabilities against the results obtained from buckling stability analysis based 
on the hyperelastic constitutive relation. FEA implementation of the failure function is discussed in Sec. (H): 
we discuss the procedure in Sec.(6.1), and illustrate its use in an example of limiting blister-type deformation 
in Sec.(6.2). Finally, we conclude in Sec.([^ by summarizing our results and their implications for analysis 
of strength-measuring nano-scale contact experiments. 


2. Framework of representation 
2.1. Kinematics 

We consider graphene as a 2 D deformable body denoted by unstressed reference configuration B. Let X 
and X denote the coordinates of a material element of graphene in undeformed and deformed configuration, 
respectively. The convection of material points under deformation is described by a smooth, injective (one- 
to-one) function x(X, t) called the motion. The non-translational part of the motion can be equivalently 
defined by the positive-definite second-order deformation gradient tensor, F = Vx(X,t). Then, the polar 
decomposition theorem provides the following factorizations of F [unann]: 

F = RU = VR, (1) 

where the orthogonal tensor R G SO 2 characterizes rigid-body rotation, and U (or V = RUR^), termed 
the right (left) Cauchy-Green tensor, characterizes shape- and area-change. The deformation of a material 
point can be kinematically factored as the product of a purely dilatational (or shape-preserving, but area¬ 
changing) deformation U^, and a purely isochoric (or shape-changing, but area-preserving) deformation U. 
Accordingly, the stretch tensor can be product-decomposed as 

u = U^U = UU^, (2) 

where 

U“ = I (3) 

and 

U = Ari (g) ri + A“^r 2 G Y 2 ] (4) 

here J = det U, A > 1 is the deviatoric stretch, I is the 2 D identity tensor, and ri = cos^ ei + sin^ 62 and 
Y 2 = — sin^ ei + cos^ 62 are the major and minor principal stretch directions, respectively. 

Graphene is a composite of two Bravais lattices shifted with respect to each other by a shift vector do- 
Upon deformation, the individual Bravais lattices follow the macroscopic deformation gradient F, i.e., the 
lattice vectors in the deformed crystal are given as = Fai and a^ = Fa 2 . This is called Gauchy-Born kine¬ 
matics. However, the shift vector, separating the two lattices, does not obey the Gauchy-Born kinematics, 
i.e., d 7 ^ Fdo- The shift vector d, in a deformed state, is determined by minimization of total energy of the 
deformed crystal. That is, under certain imposed deformations, graphene experiences sub-lattice shifts that 
lower the total energy of the crystal compared to the Gauchy-Born unrelaxed crystal. These deformation 
modes are the ones that involve a deviatoric stretch, i.e., U 7 ^ I. Examples of this class of deformations 
include the uniaxial stress/stretch and simple shear. Gonversely, because of symmetry, a purely volumetric 
deformation does not generate a sub-lattice shift. 

The sub-lattice shift s, associated with a deformation gradient F is defined as the difference between the d 
vectors connecting the two Bravais lattices in the relaxed and unrelaxed configurations. Mathematically, 

s = d - dcs; where dcB = Fdo, (5) 
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where do denotes the shift-vector in the undeformed configuration. 

The rigorous energetic and kinematic continuum description of a complex crystal, such as graphene, requires 
the definition of the strain energy density function as follows 

Ip = 'ip(F,s). (6) 


The above description implies the following variational form 


d'lp = 


^:SF + ^ 

dF ds 


Ss. 


Localization of the principle of virtual work then requires the following energy balance 


0 = 


dilj 


_ 



: Js, 


(7) 

(8) 


where is the first Piola-Kirchhoff stress (stress tensor work-conjugate to F), and f are 
unit reference volume acting on the atoms of the state defined by (F, s). The conditions 
equilibrium equations 


dF 


the forces per 
imply the 

(9) 


and 



( 10 ) 


At equilibrium, the force on each atom in the lattice is zero, and the above equation becomes 



( 11 ) 


Thus equation implicitly determines the equilibrium value of sub-lattice shift, s = Seq(F), for which 
the crystal satisfies both external and internal equilibria. Inserting this form into the strain energy function 
allows representation of the equilibrium strain energy in terms of F alone: 


^ = ^(F,Seq(F))=^(F). 


( 12 ) 


2 . 2 . Strain measure 

In this work, the logarithmic strain measure = InU is employed in the representation of the stress- 
strain constitutive response as well as of the limit surface. The spectral representation of = InU = 
In + In U is given by: 


E(o) 


- In JI In A (ri (g) ri - r2 (8) r2) 

A V._ _ ^ 



where 

ea = trE(°) = In J = ln(det U), 
gives the areal logarithmic strain e^, and 


= In U = In A (ri (g) ri - r2 (8) r2), 


(13) 

(14) 


(15) 
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denotes the deviatoric part of The orientation 0 is measured from a fixed material axis aligned with 

the zigzag direction of the graphene lattice. In the rest of this article, we will refer to 0 as the deviatoric 
stretch angle. Alternatively, we may also write the spectral representation of the logarithmic strain as 

=^maxri 0ri+£:minr2®r2. ( 16 ) 

where fmax = ln(J^/^A) and fmin = ln(J^/^/A) are the major and minor principal logarithimic strains, 
respectively. In addition, we also define a mean hydrostatic strain £ = ^(fmax + ^min) = In which 
characterizes the areal dilatation or contraction of the material. 

2 . 3 . Symmetry-constraints on the limit surface 

The proposed approach is based on the notion of a limit function — a continuous scalar-valued function 
T of the strain measure — such that all deformed configurations of the lattice that lie on the surface 

= 0, (17) 

have reached a state of incipient lattice instability. All deformed states of the crystal lying in the interior of 
the limit surface satisfy > 0, and are stable with respect to lattice perturbation of all wavelengths; 

conversely, for deformed states lying outside the surface, i.e., < 0, the lattice is unstable with respect 

to an incremental perturbation of some wavelength. 

For an anisotropic material, the limit function T — in accordance with Neumann’s principle [ 33 ] — should 
include the material symmetry group Q of the underlying lattice, i.e., 

^■(Q^E^o^Q) = V Q e a, ( 18 ) 

where Q is an orthogonal tensor denoting the symmetry operations included in the material symmetry group 
Q {Csy in the case of graphene). A scalar-valued function of a tensor agency— such as J^(E*^^^) — that remains 
invariant under a material symmetry group Q is called a ^-invariant scalar function. The representation of 
a generic ^-invariant scalar function involves using the isotropicization theorem and symmetry-invariants of 
the tensor agency with respect to the point symmetry group of the material. 


2 . 4 - Isotropicization theorem and symmetry-invariants 

The isotropicization theorem — based on the notion of a materially-embedded structure tensor H — 
allows a ^-invariant scalar function to be expressed in terms of a list of special scalar functions — J7i, ^2, 
..., Jn — which are joint isotropic functions of E^^^ and H [HISlISl]; i.e., 

^(E(°))=j^(Jl,J2,...,Jn), (19) 

where 

Ji(E( 0 ); H) = Ji(Q^E( 0 )Q; Pq(H))V Q G SO2. ( 20 ) 

Here Pq denotes the transformation of the structure tensor H under the orthogonal transformation Q. The 
functions ffi are termed symmetry-invariants since they satisfy all the constraints belonging to the material 
symmetry group of the crystal. Smith [ 5 ^ [ 53 l [ 53 ] showed that the set of mutually-independent symmetry- 
invariants serves as a complete and irreducible basis for the representation of scalar constitutive functions 
of the anisotropic material. In the following section, we explicitly derive symmetry-invariants of E^^^ for the 
structure tensor characterizing the material symmetry group of graphene. 
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2.5. Symmetry-invariants and functional basis for material symmetry group 



Figure 1: (a) Graphene lattice with orientations of the material unit vectors — x and y — and the Cartesian unit 
vectors — ei and 02 — indicated. The dashed blue lines denote the unit cell used in the ab initio calculations. The 
GGA (LDA) value of the lattice parameter is also indicated. The armchair and zigzag directions are along the ei 
and 02 axes, respectively, (b) Brillouin zone of graphene with high symmetry points and irreducible wedge indicated. 

The structure tensor characterizing the Cqv point group, the material symmetry group of graphene, is a 
sixth-order tensor given by [HSl US] 

EI = M(g)M(g)M-(M(g)N(8)N + N(g)M(8)N + N(g)N(g)M), (21) 

where M and N are two traceless second order tensors given by 

M = x(g)x-y(g)y; N = x(g)y + y(g)x. (22) 

X and y denote orthogonal material unit vectors fixed in the frame of the reference crystal such that at least 
one of them is aligned with an axis of reflection symmetry. 

The complete and irreducible set of polynomial joint invariants of and H constitute what is called an 


integrity basis, and are given by 

Ji =ea=trE(0) =lnJ, (23) 

1 : e[,°) = (In A)^, (24) 

where A : B = tr(A^B) is the scalar tensor product, and 

Ja = {le/2f = 1h[e(°\ E[,°t e[,°)] (25) 

= 1 (M:E[,°))'-3(M:E(°)) (N:E[,°))' (26) 

= (InA)^ cos66>, (27) 


where cos^ = ri.x indicates the orientation of maximum principal stretch. The first two of these invariants, 
Ca and 7 ?, are simply two isotropic invariants of alone. Thus any material anisotropy in the constitutive 
response of graphene is captured solely by the third invariant 7 ^. In a previous work [25], we showed that a 
constitutive function of graphene — such as its finite deformation hyperelastic response — is conveniently 
represented in terms of the integrity bases of the logarithmic strain with respect to the material symmetry 
group of graphene. 
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For the purpose of representation of T^ we employ an alternate set of symmetry-invariants of com¬ 
prising the mean hydrostatic strain f, the magnitude of deviatoric strain 7 ; and the symmetry-reduced 
principal stretch direction defined as 0(6>) = arccos(cos 66 >). Note that the elements of this set are not inde¬ 
pendent from the elements of the integrity basis, i.e., e^, 7 ^? 

equations ( T3p5 ), the two sets of symmetry-invariants can be shown to be related as follows 


and 7 I; and utilizing kinematical definitions of 


1. Mean hydrostatic strain f = ^ = In 


(28) 


/E(o) : E(o) 7i 

2. Magnitude of deviatoric strain 7 = \/-^- = — = | In A| >0. 


(29) 


3. Directionality 0(^) = arccos 


7 | 


= arccos(cos 6 ^). 


Thus, the failure function has the representation 


(30) 


(31) 


The set of symmetry-invariants, 7 , E and 0(^), which involve non-polynomial combinations of strain com¬ 
ponents, is called a functional basis. 


2.6. Features of the limit surface 

The limit surface — as expressed in terms of the elements of the functional basis — prescribes the critical 
value of the magnitude of deviatoric strain 7 "^ as a function of the areal strain £ and the symmetry-reduced 
principal stretch direction 0(^). Thus, the form of the limit function reads as 

^(E(o))=7<=(f,0(0))-7; 7>0, (32) 

such that if, at a given state of deformation, 7 ^(E*^^^) > 0 , then the material is stable with respect to incre¬ 
mental perturbations of all wavelengths; otherwise the instantaneous configuration of the lattice is unstable 
under some incremental perturbation. The limit surface, when graphically expressed in terms of f, 7 and 
0(^), constitutes a 3D polar surface with £ on the radial axis, 7 on the vertical axis and 0 as the angular 
coordinate. 

Boundedness of the limit surface — The proposed limit surface is bounded both on the vertical axis 
(i.e., 7 -axis) and the radial axis (i.e., f-axis), as explained in the following. 

On the vertical axis: 

• By definition, the magnitude of deviatoric strain is a positive semi-definite entity, i.e., 7 > 0 at all 
values of £. Thus, the limit surface is bounded from below by the surface 7 = 0 . 

• From above, the failure surface is bounded by the critical value of the magnitude of deviatoric strain 
beyond which the crystal is unstable, i.e., 7 < 7 ^(f, 

On the radial axis: 

• Because of graphene’s limiting thickness and extreme bending compliance, essentially any attempt 
to decrease the area of a sufficiently large, flat graphene sheet would result in out-of-plane buckling, 
graphene’s small, though non-zero bending stiffness notwithstanding; for this reason, we consider that 
the domain of the surface is bounded from below: f > 0 , providing a graphical representation using 
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f as a ‘radial’ coordinate. Moreover, if £ is held at 0, imposition of a non-zero 7 creates a negative 
principal stress, again leading to buckling in a sufficiently large, flat graphene sheet; hence we take 
7^(5 = 0, 0(6>)) = 0 for all 0. (These limits follow from the quadratic phonon dispersion relations of 
ZA phonons at T in a periodically-extended planar crystal.) 

• Upon a continued equi-biaxial deformation, a state is reached, £ = fmax when graphene fails without 
any shear strain, i.e., in pure dilatation. Therefore, the radial domain of the surface is bounded from 
above as well: £ < fmax and 7 ^(f = fmax, 0 (^)) =0 for all 0. 


3. Form of the limit surface in terms of the functional basis 

Employing the representation ideas outlined in the previous section, we now systematically determine 
the functional form of the limit function T = 7 ‘^(f, 0 (^)) — 7 for graphene in terms of the three symmetry- 
invariants. For this purpose, we first obtain the set of data containing the stability-limiting values 7 "^ 
over a range of £ values, calculated from lattice-dynamical stability analysis, for a set of representative 
homogeneous deformation paths. 


S.l. Lattice dynamical stability analysis 

The atoms in a crystalline material remain in a state of oscillatory motion about their equilibrium sites, 
constituting what is called the random lattice vibrations. Phonons are the normal modes of the lattice 
vibrations. Fundamentally, each phonon is a collective oscillation of atoms characterized by a wave-vector 
q and a oscillation frequency where n labels the branch index. The relation between and q is called 
the dispersion relation of the crystal. The phonon-dispersion relation of a crystal contains a total of SA/" 
branches where M is the number of atoms in the unit cell {M = 2 for graphene). Three (3) of these branches 
are called acoustic and the remaining 3(A/' — 1 ) are called optical. 

Phonons constitute a useful means of assessing the mechanical stability of a crystal lattice. Mechanical 
stability of a crystal requires that the underlying lattice remains stable against spatial perturbations of all 
wavelengths, and this is ensured if the following two conditions hold: 

• First, all the acoustic branches should have a positive slope at q = 0 , i.e.. 


>0 VnG 1,...,3. (33) 

q=0 

If at some point along the loading path, this criterion ceases to hold then it is an indication of an 
incipient instability, called a long wavelength instability, in the material. In monolayer graphene, 
a long wavelength instability may correspond to either an out-of-plane instability (buckling) or an 
elastic instability Buckling is a mode of structural instability wherein graphene explores the 

third spatial dimension via deformation under compression. Unlike a material instability, buckling 
does not involve material fracturing or undergoing plastic deformation. Further, since graphene is one 
atomic layer thin, a buckling instability manifests itself as a long-wavelength instability in out-of-plane 
acoustic phonon (also called ZA branch) mode of graphene. 

On the other hand, an elastic instability is a mode of material failure and manifests itself as a long- 
wavelength instability in the in-plane acoustic phonon (LA or TA branch) modes of graphene. An 
elastic instability in the material can also be captured by acoustic tensor analysis; therefore the onset 
of elastic instabilities in material can be parameterized with strain via the acoustic tensor route also. 

• Second, all phonon frequencies for q 7 ^ 0 should be real and non-zero, i.e., 

cc;^(q 7^ 0) > 0 V n G 1 ,..., 3A/'. (34) 


duJn 

(i|q| 


If at some point along the loading path, this criterion ceases to hold then it is an indication of an 
incipient soft mode instability in the material. Unlike an elastic instability, a soft mode instability does 
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not have an equivalent continuum criterion, and therefore a mathematical framework parametrizing 
the onset of such instabilities in terms of stress or strain. 

Lattice-stability analysis based on the above criteria is most effectively assessed by phonons, which, below the 
Debye temperature (Td = 1000LC for graphene), constitute a complete normal basis for lattice vibrations. 
A deformation-induced instability of any kind — macroscopic or microscopic — is directly visible in the 
dispersion of phonons. 


3.2. Basics of density functional perturbation theory 

DFPT employs density functional theory (DFT) in conjunction with linear response perturbation theory 
to calculate the lattice dynamical properties of a crystal (see references [2j [T31IH] ). In density functional 
theory (DFT) [^[26] . the ground state energy of a system of interacting electrons moving under the influence 
of an external potential is obtained by minimization of a universal functional, J^[n(r)], of electronic density 
n(r). The electronic density no(r) that minimizes the energy functional also happens to be the true ground 
state electronic density of the system. 

Kohn and Sham m showed that the original system of interacting electrons can be replaced by an auxiliary 
system of non-interacting electrons moving under the influence of an effective potential such that the auxiliary 
system at the ground state possesses the same electronic density as the original system. The motion of 
electrons in such an auxiliary system are described by a set of one-electron equations (also called Kohn- 
Sham equations): 

V’n(r) = enV’n(r), (35) 

where ipri is a Kohn-Sham orbital and is the corresponding eigenvalue. The electronic density is calculated 
from the Kohn-Sham orbitals using the relation: 

= '^\'4’n{r)\'^g{eF - ei), (36) 

n 




where ^(ei? — e^) is the occupation function, and ep is the Fermi energy. The effective potential KEff(i*) is a 
functional of electronic density, and is given by 


I4:ff(r) 


yiE(r) + y \r-r'd ^' 


(37) 


where the first term on the right, IdE(r), is the potential due to interaction between electrons and ionic 
cores, the second term is the potential due to Coloumb interaction between electrons, and the last term, 
'^ccc(i*) = ^^ccc/^^(r), is the functional derivative of the exchange-correlation energy functional Sxc- Once an 
explicit form for £xc has been determined, the set of equations (35 - 37) can be solved self-consistently to 


determine the ground state energy and electronic density of the system. 

Calculation of phonons requires knowledge of interatomic force constants, defined as the derivatives of the 
ground state total energy of the system E with respect to ionic coordinates R = {Ri,R 2 , ...,RAr}- Using 
Feynman-Heilman theorem, we evaluate such derivatives as: 


a2F;(R) 

5R/9Rj 


/ 


dn{r) 

aRj 


«o(r) 


51^iE(r) 

aR/ 


dr + 



d^Vmjr) 

dUidHj 


dr + 


a2Uii(R) 
dKidKj ’ 


(38) 


where Un is the energy corresponding to Coloumbic interaction between ionic cores. Thus, the evaluation 
of the interatomic force constants requires ground state electron density no(r) as well as the first-order 
correction [dn{T) 

The DFPT method calculates the first-order correction in electronic density from lattice’s response to a 
set of monochromatic lattice perturbations, each characterized by a q- vector. Provided the amplitude of 
such perturbations are small enough, the application of perturbation theory allows to recover a set of self- 
consistent relations for the first-order corrections in electronic density and wavefunctions. Let and 
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be the unperturbed wave function and eigenenergy of an electron with wave-vector k and band-index u, 
respectively. First-order correction in the wave function then satisfies the following equation: 


- V^ + VEff 
2m 




(39) 


where is the projector operator over the manifold of states of wave-vector k -t- q, and AVEff(r) is the 

associated perturbation in the effective potential due to perturbation, given by: 


AyEff(r) = J 


^n{r') 3 


r - V 


+ An(r) 


dv^c(r) 


dn 


■ AhiE(r). 


no(r) 


Perturbation in the electronic density is obtained as: 


An = - 


(40) 


(41) 


V being the volume of the system. Upon solving the set of equations ( 39pT ) self-consistently, we obtain the 
first-order corrections in wave function and electronic density of the system. Knowledge of such corrections 
allows direct access to the dynamical matrix of the system, defined as the Fourier transform of the interatomic 
force constant matrix, i.e., 


D)/aJ6(q) 


1 




y^C/Qjb(R)e 


— iq.R. 


ClaJb 


R 


d^E{K) 

dKiadK 


Jb 


(42) 


where M/(j) denotes the mass of atom. Phonon frequency Cc;(q) and the associated eigenmode Uq 

are then determined by solving the eigenvalue equation: 


CJ^(q)Uq = D(q)Uq. 


(43) 


3.3. Details of DFPT calculations 

The exchange correlation energy of electrons is treated with the Local Density Approximation (LDA) 
of Perdew and Wang (ESI). The interaction between ionic cores and valence electrons is represented by an 
ultrasoft pseudopotential m- Kohn-Sham wave functions are represented using a plane-wave basis with 
an energy cutoff of 30 Ry and a charge density cutoff of 300 Ry. Integration over the irreducible Brillouin 
zone (IBZ) is performed with a uniform 30 x 30 x 1 mesh of /c-points, and occupation numbers are smeared 
using the Marzari-Vanderbilt cold smearing scheme with broadening of 0.03 Ry. Errors in the Cauchy 
stresses and total energy due to basis-set size, smearing parameter, and /c-points are converged to less than 
0.034 N/m and 0.01 Ry, respectively. 

The phonon dispersion relations — of the undeformed and deformed graphene — are computed via linear 
response calculations as implemented in the density functional perturbation theory (DFPT) [2l[T3l. The 
dynamical matrix is calculated on an 8 x 8 x 1 uniform grid of q-points in the IBZ — which is then fast- 
Fourier-transformed to calculate the interatomic force constants (IFC), corrected by the acoustic sum rule 
to ensure that Cc;(q = 0) = 0 for all the acoustic branches. The IFC’s are then used to interpolate phonon 
frequencies over a dense set of q-points in the IBZ. Both energies and phonon dispersion relations in this 
work are performed on a two-atom primitive unit cell of graphene shown in Fig. 0. All calculations are 
fully relaxed in terms of shift vector to ensure conditions of zero atomic force. 

The LA and TA phonon branches, in the neighborhood of q = 0 , obey linear dispersion: ujla = claQ 
and ujta = ctaQ] where g = |q| and cla and cta are longitudinal and transverse acoustic wave-velocities, 
respectively. As a simple verification of the DFPT phonon calculations. Tab. 0 compares longitudinal 
and transverse acoustic wave-velocities obtained from the phonon dispersion with experimentally-measured 
values; also, since the G band in the Raman spectra of graphene is associated with the doubly-degenerate 
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LO - TO vibrational modes at T, we also compare calculated LO/TO frequencies at T with the G-band 
frequency measured in Raman spectroscopy. 



CL^(km/s) 

CTA(km/s) Raman G-peak 

Measured value [58] 

21.3 

13.6 1582 cm-i 

Continuum model (LDA) [25] 

21.80 

13.78 — 


V E2dIP2D 

V P2dIP2D 


DFPT Phonons in this work (LDA) 


21.15 


duJLA 


dq 


q=0 


13.50 


duJTA 


dq 


q=0 


1540 cm ^ 


Table 1: In-plane longitudinal and transverse wave velocities calculated from phonon dispersion closely 
agree with measured values as well as with values calculated from the continuum model [25]. The LO/TO 
frequencies at F are also in good agreement with the G-peak of the Raman spectra of graphene. Note: 
P 2 d is the mass per reference area of graphene which is calculated from the atomic mass of carbon and 
the undeformed lattice constant of graphene; E 2 D denotes the elastic modulus in uniaxial strain; and p 2 D 
denotes the shear modulus. 

DFPT allows calculation of phonon dispersion relations with high accuracy at any arbitrarily-deformed 
state of the graphene lattice, whereas the errors in phonon frequencies based on empirical potentials can be 
forbiddingly large, rendering the predictions somewhat unreliable. For example, the sound velocities from 
an empirical many-body potential in the unstrained state of graphene differed from measured values [55] 
by nearly 50% [12]; in contrast, the longitudinal and transverse elastic wave-velocities obtained from DFPT 
phonon calculations differ by less than 1% from measured values, as shown in Tab. 0 . 


3.4- Sampling scheme 

For the purpose of sampling deformed states of the lattice, we use the deviatoric stretch angle 0 to 
designate what we term a ‘deformation path’; i.e., increasing 7 at a fixed 0 is referred to as moving along the 
deformation path prescribed by To determine the point of incipient lattice instability along a deformation 
path: 


• The lattice is first deformed via a pure equi-biaxial strain, i.e., = El , followed by an isochoric 

shape-changing stretch of the form = 7 (riG)ri — r 2 ( 8 )r 2 ) where ri = cos 6 >ei+sin 0 e 2 and r 2 T ri; 
0 is kept fixed. Thus, total strain in the final configuration is E^^^ = E^^^ + Eq^^ , as schematically 
shown in Fig. 0 a). 


For each deformation pair {5, 7 } along the sampled deformation path 6 >, we check the phonons to 
determine if the conditions for acoustic stability (given by equation (33)) and short-wave stability 
(given by equation @) hold. If at some critical magnitude of deviatoric strain, 7 = 7 "^, either 
of the two conditions ceases to hold then the lattice is considered as incipiently unstable, and the 
corresponding triplet {5, 7 ^, 0} corresponds to a point on the instability surface. 


Our sampling set includes values of £ ranging from the undeformed state £ = 0 to the critical equi- 
biaxial deformation, £c = 0.144, at which the lattice reaches a soft mode instability in the absence of 
deviatoric deformation. For each sampled value of 5, and along each sampled deviatoric stretch angle 6 >, 
the superimposed deviatoric strain is varied such that its magnitude 7 varies over the range 0 < 7 < 7 "^, 
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the upper limit 7 ^ being the critical value of the stretch at the given value of £ and 0. Owing to the 
symmetry of graphene, the deviatoric stretch angle 0 needs to be sampled only over the range 0 < 6 > < 7 r/ 6 . 
Further, we consider only N = 2 deformation paths in the BZ, corresponding to the zigzag direction ^ = 0, 
and the armchair direction 0 = tt/G. Later, we will show that only two sampling directions suffice for an 
accurate angular interpolation of the failure function. 

3.5. Representation and interpretation of the data 

From the lattice-dynamical stability analysis, we obtain the limiting values, each in the form a triplet 
7 ^, 6 >}, representing the critical value of the deviatoric strain magnitude — 7 ^ — as a function of mean 
hydrostatic strain £ for both the sampled deformation paths: 0 = 0 (shown in Fig. §b)) and 0 = 7T16 
(shown in Fig. §c)). 



Figure 2 : (a) Notation scheme illustrating deformation of the graphene lattice: the undeformed lattice (faint black) 
is first subjected to a uniform dilatation Ui = (resulting intermediate configuration shown in faint blue), and 

then to an isochoric shape-changing deformation U2 = Ari (8)ri + A“^r2(8)r2 (resulting final configuration shown in 
dark red), (b) Cartesian plot of 7^ as a function of £ for the sampled orientation ^ = 0 , which corresponds to the 
zigzag direction, (c) Cartesian plot of 7^ as a function of £ for the sampled orientation 0 = tt/G, which corresponds 
to the armchair direction. 


Shown in Fig. (§ are the phonon dispersion relations of graphene calculated with DFPT at certain critical 
deformed states, indicated as ‘1’, ‘2’ and ‘3’, along the sampled deformation path ^ = 0; these particular 
critical states of deformation are those labeled on Fig. §b). Depending upon the level of mean hydrostatic 
strain, a critical phonon state might be associated with buckling, as in ‘ 1 ’, for which the long-wave ZA 
phonon slope at F vanishes; with an elastic (material) instability, as in ‘ 2 ’, for which the long-wave LA slope 
at F vanishes; or with a short-wave material instability as in ‘3’, where the frequency of the TA phonon van¬ 
ishes at K. Thus, at different states of deformation, different kinds of mechanisms govern lattice instability 
in a graphene sheet. The proposed limit surface is a smoothed envelope of all the possible lattice-instabilities, 
as we will elaborate further in Sec.(3.7). 
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Figure 3: Phonon dispersions of graphene at certain critical states of deformation, marked as ‘1’, ‘ 2 ’ and 
‘3’ in Fig. it), along the armchair deformation path (^ = 0). Also shown are the deformed Brillouin zone 
(green polygon), and the irreducible wedge (red polygon). The sampled direction in each case is shown by 
an arrow. The labels ‘ZA’, ‘LA’ and ‘TA’ show the branch that goes unstable in the three cases. 


3.6. From discrete dataset to continuum lattice-stability limit function via interpolation 

Construction of a continuum-level limit function involves interpolating the discrete dataset of the pre¬ 
ceding section in a manner that is consistent with the material symmetry {CQy) of graphene. This task is 
accomplished by interpolating the critical magnitude of deviatoric strain 7 *^ in terms of two of the symmetry- 
invariants of the strain tensor: £ and Q{0). The interpolation is performed in two steps: first, interpolation 
in terms of f, called radial interpolation, and second, interpolation in terms of 0 (^), called angular inter¬ 
polation. 

3.6.1. Radial interpolation (interpolation in terms of £) 

Along each sampled deformation path 0 = 0^, we express 7 ^ by a continuous function lZn{£) — called 
a radial interpolation function. Each such function is obtained by fitting the sampled datapoints (plotted 
in Fig. (|^a)), which are in the form of a doublet , 7 *^} and denote the limiting strains along a sampled 
deformation path, to a polynomial function in £ of the form: 

7"(^,0 = 0n)=7^n(^/fc) 

= + Pl{£/£cf + + ... + P^{£l£c)^, (44) 

where fc = In Jmax = 0.144 is the mean hydrostatic strain at which an incipient lattice instability emerges 
without any deviatoric strain. The radial interpolation functions thus obtained, corresponding to the two 
sampled deformation paths, are shown in Fig. @b). 
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Figure 4 : (a) Sampled valued of 7^ along two deformation paths: ^ = 0 and ^ = tt/G (straight-line segments connect 
data points), (b) 7 ^^-order polynomial interpolation functions for 7^ along the armchair and zigzag directions — 
obtained by htting functions described by equation 


44 


3.6.2. Angular interpolation 

The angular interpolation scheme involves approximating the failure surface 7 ‘^(f, 0) by a weighted sum 
of the radial interpolation functions — with the weight functions being the angular shape functions along the 
sampled deformation paths. The following representation for the angular shape function /Cn(0) conveniently 
satisfies the requirements that failure function is smooth, satisfies all the point-group symmetries of the 
lattice, and possesses 27r periodicity with respect to 0 : 


A^n( 0 ) = <^n + cos( 0 ) + a‘^ cos( 20 ) + ...; 


(45) 


where the number of terms in the expansion is given by the number of deformation paths (N) being sampled, 
and the corresponding coefficients are determined using the following condition: 


JCmiQn) 


0; 0n 7^ 0m 
I5 0n — 0m* 


(46) 


Since the angular variation of the limit function is quite smooth, sampling along only two deformation paths, 
i.e., N = 2^ suffices for accurate interpolation of the limit surface (see Fig. in Sec. ([^). The corresponding 
two shape-functions are given as: 


^i(©) = ^(i + cos(e)), 

(47) 

/C2(e) = l(i-cos(e)). 

(48) 


Following the above interpolation scheme, 7 ^ is generically expressed as: 


n=N 


Y{£,o{ 9)) = iiu{£/£c)iCn{em, 


(49) 


n=l 


where = 2 is the total number of deformation paths. Substituting from equation (44) and equation (45) 
into equation (49), we obtain the generic form of the lattice stability function as: 

n=N M 


^%£,eie))= 7mn(74)™cos(n-1)0(0); 


(50) 


n=l m=l 
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The various coefficients appearing in the expression of the failure function have been tabulated in the 

Tab. (|§. 


Jmn m ->■ 

n 4 . 

£l£c 



77c)" 




1 

0.21561 

- 0.0403958 

- 0.327806 

7.51334 

- 24.8089 

28.3255 

-10.8837 

COS 0 

0.0251542 

- 0.583471 

4.55324 

-16.0798 

26.6759 

- 20.8193 

6.22982 










Table 2: Coefficients ^mn in the expression of the failure surface given by equation (50); the basis function 
for each coefficient 7 ^^ is the product of the column of the top-most row with the row of the 
left-most column. 


The 0 < ^ < 7r/2 quadrant of the instability surface resulting from the interpolation scheme is shown in 
Fig. (|5}b). 



Figure 5: (a) The radial interpolation of the lattice-stability limits along the two sampled directions, ^ = 0 and 
0 = TV 12. (b) A 90^ quadrant cut out from the limit surface. From below, the limit surface is bounded by the zero 
surface shown in green; while from above the failure surface is bounded by the critical surface 7 = 7 ^(^, C(^)), shown 
in red. The directions ^ = 0 and Q — j2 correspond to the zigzag and armchair directions, respectively. 


S. 7. Limit surface as a smoothed envelope of various possible instabilities 

Different kinds of strength-limiting mechanisms lead to instability of the graphene lattice under different 
modes of deformation. The limit surface is a smoothed representation of the envelope of all possible lattice 
instabilities: long wavelength as well as short wavelength, and structural as well as material failures. 

Graphene, being a strictly two-D structure, has an exceedingly small bending rigidity. Therefore, when 
subjected to an in-plane deformation, a graphene sheet can be brought to a state of incipient instability via 
one of two mechanisms: 

• a material instability in tension, which arises when the major principal strain fmax (maximum strain 
over all material directions) reaches a limiting value, while the minor principal strain Emin (minimum 
strain over all material directions) has not reached the critical value in compression so that both the 
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Cauchy principal stress components are non-negative, i.e., iSmax > <5niin > 0. This material instability 
can be either of the two types: long-wave (also called elastic) or a short-wave (also called soft mode) 
or 

• a structural instability via out-of-plane buckling, which arises when the minor principal strain fmin 
(minimum strain over all material directions) in the material reaches the critical value in compression 
such that one of the principal stress components tSmax or ^^min becomes incipiently compressive (< 0 ), 
while the major principal strain has still not reached the limiting tensile strain. 

The two principal strains, in terms of £ and 7, are given by fmax = f + 7 and fmin = f — 7 . Thus, at 
a fixed equi-biaxial strain f, with a superposed, increasing deviatoric strain 7, the major principal strain 
monotonically increases while the minor principal strain monotonically decreases. Now, whether a deviatoric 
strain 7 superposed on an equi-biaxial strain £ results in a material instability or a structural instability 
depends on the magnitude of 5, as schematically illustrated in Fig. For example, when £ is small, 
the minor value £min reaches the critical value in compression, which is ~ —z^ 2 i:>^max, before the major 
principal value fmax reaches the limiting tensile strain. Thus, at small equi-biaxial strain, the magnitude 
of the deviatoric strain 7 is essentially limited by a buckling instability. This is also shown in the phonon 
dispersion relation of Fig. § 1 ). Interestingly, when f = 0, the membrane is unstable in buckling as soon 
as any deviatoric strain is applied; thus the limiting value of 7 ^|^^o = 0 V 6 >. 

On the other hand, when £ is large, fmax reaches the limiting tensile strain before £min reaches the criti¬ 
cal value in compression, and the failure is invariably a material failure (fracture). Further, as £ approaches 
closer to the critical value £c, the material starts to fail essentially in dilatation via a short-wave instability, 
as shown in Fig. (|^3). Thus, ^ 0 V 



Figure 6: Schematic illustration showing the lattice-stability plot as a smoothed representation of the envelope of 
the various possible lattice instabilities. 


4. Lattice-stability limit surface in stress space 

Often, it is useful to describe the lattice-stability limits in terms of stress. In the present case, although 
the limiting conditions are explicitly obtained in terms of the invariants of strain, the relationship can be 
numerically mapped to its counterpart in the stress space via a constitutive function appropriately describing 
the stress-strain response in graphene over the entire range of deformation up to the elastic stability limit. In 
a previous work [25] , we derived a large-deformation hyperelastic constitutive response function for graphene 
based on ab initio calculations — which we briefly describe in the following. 
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The constitutive response function is derived within the framework of nonlinear hyperelasticity and invariant- 
based representation theory. Our formulation employs the symmetry-invariants of the log strain measure 
— ea , 7 ?, and 7 I — to express the dependence of the free energy per unit reference area ip on the 
state of strain i.e., pj = '0(6^, 7 ^^, 7 I) (see Kumar & Parks [28] for details). The strain energy density 
function (e^, 7 ^, 7 !) accepts the additive decomposition: p; = 7 ^^, 7 I). The term 

— the energetic response under pure dilation — is described by a function based on the universal binding 
energy relation gaiizi: 

'Ip^'\ea) = ■^o [1 - (1 + aea) exp(-Q:ea)]. (51) 

The term — the energetic response under an isochoric deformation — is given by a simple additive 
form: 

= \K<^ahf + g??(ea)7e> (52) 

where = /j^q — is a second order elastic constant, and r]{ea) = Vo — is a third-order elastic 

constant. The coefficients are determined from fits to ab initio energies for a set of deformed configurations 
(see Tab.[^). Differentiation of pj with respect to E*^^^ gives the work-conjugate stress tensor, i.e., = 


t/)o(J/m2) 

a /io(N/in) 

Mi(N/m) 

13 

%(N/in) 

7?i(N/m) 

LDA 116.43 

1.38 172.18 

27.03 

5.16 

93.17 

4408.76 


Table 3: Coefficients in the expression of strain energy function '0(ea, 7 ^, 70 ). 


dpj/dE^^\ 


7 ^ 7?) = 

5V;^^^(ea,7^^7|) U , {ea^hle)-^{o) , {ea.-fhle) ^ 

v^* m) " m) ■="’ 

= |v'oti^<^aexp(-aea) + ^p'(<^a)7i + ^V{ea)7l j I + 2p(ea)Ej,^’ + ^j){ea)Sg^o), (53) 

and the Cauchy stress a is given as (see Ogden [33]) 


<T = = 2 jF 


dC 


rp(0)pT 


(54) 


where C = denotes the right Cauchy-Green stretch tensor. Employing equation (53) and equation (54), 
we numerically map the contours in the strain plane to stress-plane. Once the Cauchy stress has been 
determined, we can obtain the principal stresses from the spectral representation of a: 


(T — <5maxSl C Si “h <SniinS2 C S 2 , (55) 

where tSmax^ and iSmin are the principal Cauchy stress components; and Si = cos^ei + sin0e2 and S 2 = 
— sin^ei + cos 0 02 are the corresponding principal directions. The orientation p is measured from a fixed 
material axis aligned with the zigzag direction of the graphene lattice. 

Employing the above transformation, we obtain the lattice-stability plot in stress space as a functional 
relation between the components of the functional bases comprising the maximum shear stress Tmax = 
(5max — *5min)/2, thc mean hydrostatic Cauchy stress a = (<Smax + 5min)/2, and the symmetry-reduced 
deviatoric stress direction 4> = arccos(cos(60)). The limit function— in terms of Tmax, d, and 4>(0) — can 
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be expressed in a form similar to the limit function in strain space, i.e., 

'r(o') = 'rmax(^> ^(0)) - 'Tmax = 0, T^ax > 0; (56) 

such that if, at a given state of deformation, T{(t) > 0 , then the material is stable with respect to incremental 
perturbations of all wavelengths; otherwise the instantaneous configuration of the lattice is unstable under 
some incremental perturbation. 


The surfaces described by equation (56) constitute the lattice-stability limit surface in the space of symmetry- 
invariants of the Cauchy stress tensor, as shown in Fig. 0 b). 


(a) 



0.00 


30 


(b) 


Constitutive 

relation 



0 



Figure 7: Numerical mapping of the instability surface from strain space (shown in (a)) to stress space (shown in 
(b)) via the hyperelastic constitutive relation of equation([^. The axes in stress space are in units of N/m. The 
directions 0 = 0 and 0 = 7r/2 correspond to the zigzag and armchair directions, respectively. 


5. Validation 


The limit surface — given by equation (50) — contains information of all possible lattice-instabilities 
that can be triggered by a homogeneous deformation. Specifically, once this failure surface is determined, 
we can assess the stability of graphene subject to a biaxial state of strain — 8^^ ri (g) ri + £22 1*2 C) r 2 — 
for any biaxiality (the Eh to £22 ratio) and any deviatoric stretch angle 0. For such deformations, the 
principal strains are directly given by the strain components £^1 and £22 as fmax = rnax[ff]^,£^ 22 ] ^rid 
fmin = rnin[ff]^,^ 22 ]* Substituting a particular value of 0 in the limit surface, we can obtain the 2D radial 
section of the limit surface corresponding to that particular value of 0. For example. Fig. shows the 2D 
radial sections of the limit surface corresponding to 6> = 0 and 0 = 7r/6, denoted by 7^ 1^=0 —7 = 0 and 
— -)/ = 0, respectively. These contours characterize the stability of a graphene sheet subjected to a 
homogeneous biaxial state of strain referred to the set of axes defined by ^ = 0 and 0 = 7r/6, i.e., ri — Y 2 
are aligned with ei — 62 (see Fig. (§a)). 

On this contour, we have indicated points corresponding to certain special homogenous deformation cases 
such as uniaxial strain and uniaxial stress along the armchair and zigzag directions : 


1. The point indicated as given by the intersection of the line 7 = f with the contour 7 ^ 1^=0 — 7 = 0 , 
denotes the uniaxial strain along the zigzag direction. 

2. The point indicated as T^, given by the intersection of the line 7 = 5 with the contour — 7 = 0, 

denotes the uniaxial strain along the armchair direction. 

3. The point indicated as Aa, given by the intersection of the line 7 = ^ ^ ^ with the contour 

— 7 = 0 , denotes elastic failure under uniaxial tension in the armchair direction. 
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4. Similarly, the point marked kz, given by the intersection of the line 7 = ^ ^ with the contour 

—7 = 0 , denotes the uniaxial tension in the zigzag direction. 

5. The point indicated by ■, given by the intersection of the line 7 = 0 with a contour 7'^|0 — 7 = 0 (V 0 ), 
denotes the limiting strain under equi-biaxial strain. 

Following the above-mentioned procedure, it is possible to obtain the limiting contour corresponding to 
any directionality 6>, and deduce the corresponding stability-limiting values of stress/strain in uniaxial 
strain/stress. 

K (a) ,, (b) 




Figure 8: (a) Cross-section of the lattice-stability limit surface corresponding to ^ = 0. Critical points correspond to 
special homogeneous deformed state such as uniaxial stress and uniaxial strain in the zigzag direction and equi-biaxial 
strain are indicated, (b) Cross-section of the lattice-stability limit surface corresponding to ^ = 7r/6. Critical points 
corresponding to special homogeneously deformed states such as uniaxial stress and uniaxial strain in the armchair 
direction and equi-biaxial strain are indicated. 


5.1. Comparison between predieted response and values ealeulated from DFPT 

First, we assess the predictive capability of the failure function from a numerical accuracy point of view. 
For this purpose, we compute the stability-limiting critical strains —for biaxial deformations referred to 
sets of axes that are not included in the sampling set — and compare with directly-calculated values from 
a phonons-based stability analysis. As shown in Fig. the agreement between the predicted stability 
limits from the failure model and the calculated values is very close. This, in particular, highlights the 
efficacy of the symmetry-invariants based representation in constitutive modeling, since the representation 
of the failure function utilized the explicit calculations only along the two symmetry-directions of graphene. 
Previously also, in hyperelastic constitutive modeling of graphene [28], we showed that the representation in 
terms of symmetry-invariants not only reduces the number of constants in the model; but it also elucidates 
the underlying hyperelastic softening behavior of the graphene lattice. 
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Figure 9: Shown are the stability limits for deviatoric stretch angles (a) 0 = tv/30 (b) 0 = tv/13 (c) ^ = tt/IO and 
(d) 0 = 27r/15. These directions were not included in the interpolation procedure. The interpolated stability limits 
(shown by solid red lines) by the model are in close agreement with the directly-calculated values (shown by black 
dotted circles) obtained by lattice-dynamical stability analysis. 


5.2. Comparison with acoustic tensor (K-) analysis — elastic (material) instabilities in tension 


In a previous work [28], we showed that both the 2D area modulus K{ea) 


dticT 12 
dCa 


In A=0 


and the 


2D shear modulus /i(ea) of graphene experience a progressive softening with increasing areal strain. The 
progressive softening of the elastic moduli is often responsible for long-wavelength material instabilities of 
the underlying lattice. The acoustic tensor derived from the constitutive response can be used to determine 
the macroscopic instabilities of the graphene lattice, arising from the hyperelastic softening. If, for some 
pair of unit vectors m and n. 


A(m, n) = (m (g) n) : A : (m (g) n) < 0, 


(57) 


then the material is elastically unstable; here A = d‘^^ldF‘^ is the acoustic tensor (see [ 28 ] for detailed 
evaluation of A). Employing the acoustic tensor-based analysis, we determine the stability-limiting critical 
strains for varying biaxiality in a biaxial deformation with stretch axes ri — r 2 aligned along zigzag-armchair 
(ei — 62 ) set of axes, i.e., 

= SiiBi (g) ei + ^22^2 (8) 62. ( 58 ) 


We plot the resulting limiting values of £22 —obtained from acoustic tensor analysis — and superpose 
them on the phonon-instability contour as shown in Fig. (10). The agreement is close, except in the 
neighborhood of the equi-biaxial deformation. For near equi-biaxial deformat ions (5^;^ ~ ^ 22)5 instability 
of the graphene lattice is not of elastic nature but corresponds to microscopic soft-phonon modes (see 
Yevick & Marianetti [36]), which can not be detected via acoustic-tensor analysis. Following a simple 
algebraic procedure, it can also be shown that in the long wavelength limit, the phonon-based criterion and 
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the acoustic tensor based condition become equivalent, and therefore in the case of elastic instabilities, the 
agreement between results of the two analyses is close. 

5.3. Comparison with buckling analysis — structural instabilities in compression 

The 2D acoustic tensor— based on the in-plane constitutive response— captures material instabilities 
of macroscopic (elastic) nature only. However, in addition to material instabilities, structural instabilities of 
non-material nature may also arise. These instabilities are the buckling instabilities ensued by compressive 
stresses (precisely speaking, when at least one of the eigenvalues of a becomes less than zero). The buckling 
modes on the instability contour are identified as the state of strains at which one of the principal stress 
vanishes. For example, in biaxial deformations with stretch axes ri — T 2 aligned along zigzag-armchair set 
of axes, the buckling stability-limits are given by strain trajectories (shown in Blue in Fig. ([1Q|)) calculated 
from the constitutive model for uniaxial tensile stress along these two special directions. 

It is interesting to see that, in the small strain limit, the two buckling stability-limiting contours are 
related to Poisson contraction along the zigzag and armchair directions. Taking advantage of isotropy at 
small strains, linearization of the constitutive relation in terms of principal stresses and strains gives 

= E2d{S^i + V2dS^22)-, (59) 

and 

^22 ~ ^2d{£22 ^2d£ii)' (60) 

Here the asymptotic value of the in-plane Poisson ratio emerging from our constitutive model is 1^20 = 0.205 
m- The asymptotic principal strain trajectories corresponding to the condition of buckling instability are 
obtained as follows. 

Buckling in armchair direction: 


cO _ cO cO ^ cO 

^22 ~ ^11 > ^22- 


(61) 


Thus, we obtain the following relation between 7 and £: 


^ /cO cO ^ ^ 2 d ) ^ 

■r = 5(£n - £») = 


Buckling in armchair direction: 


cO _ cO cO cO 

Oil ~ ~^2D 022', ^22 ^ ^ 11 * 


(62) 


(63) 


and again, we get the same relation between 7 and £ as in equation (62) for buckling instability. The buck¬ 


ling stability-limiting trajectory on the 7 — f plane in the asymptotic limit corresponds to a straight line 
emanating from origin with slope (1 + y 2 D)l^ — '^ 2 d\ and as noted previously E 2 D is the elastic modulus 
of graphene in uniaxial strain. 
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Figure 10: (a) The 2 D radial section of the limit surface corresponding to ^ = 0 (zigzag direction), i.e., 7 = 
7 ^ 1 0 = 0 , superposed onto the elastic-stability limit curve obtained from A— analysis and the buckling-instability curve 
obtained from Poisson analysis, (b) The 2D radial section of the limit surface corresponding to 0 = tt/G (armchair 
direction), i.e .,7 = superposed onto the elastic-stability limit curve. Intersection of Poisson analysis and the 

acoustic instability analysis marks the uniaxial tension failure. Note that in large, nearly equi-biaxial deformations, 
the phonon-instability precedes the elastic instability (conforming with Yevick & Marianetti’s [36] calculations), 
indicating that in such deformations, lattice instability is due to a short-wavelength instability. 


6. Implementation of the limit surface in a continuum FEA scheme 


6.1. Outline of the methodology 

In the following, we outline a systematic procedure for implementation of the proposed limit function 
in a finite element analysis in order to provide an on-the-fiy assessment of incipient material failure during 
monotonic loading. 


First, using the continuum-level FEA based on the nonlinear constitutive model of equation (53) and 
equation (54) — we obtain the macroscopic deformation field associated with the given boundary-value 
problem at the current load/displacement increment. Note that the stress/strain field in general is 
inhomogeneous; however, if the spatial variation of elastic fields varies sufficiently slowly, each material 
element in a sufficiently fine mesh can be fairly assumed to be in the state of homogeneous deformation, 
with the local stress/strain determined from the macroscopic deformation field. 


• Based on the local stress/strain value, we make on-the-fiy assessment of the failure function at the 
material elements in the mesh for the current load/displacement step. For each increment, the material 
elements with the minimum values of the failure function are identified as the ‘weakened-spot’. 
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• The value of the critical macroscopic load/displacement is determined by identifying the load/displacement 
increment along the loading path at which the failure function at the weakened-spot first vanishes. 
The weakened-spot is then declared as a Tailed region’, and the corresponding load/displacement can 
be considered as failure load/displacement. 

• We determine value of the critical stress/strain by locating the point of intersection of the stress/strain 
trajectory traced by the Tailed’ material elements during the course of loading with the failure surface. 

• Finally, we check if the failure surface precedes or coincides with the acoustic-instability surface at the 
point of intersection. This allows us to determine whether the instability is a soft mode instability or 
an elastic instability, respectively. 

Remark 1 : Since the phonons are defined with respect to an ideal infinite crystal without any surface 
or boundary, the application of phonons to stability detection in real boundary-value problems requires a 
little care. For this we need an volume/area element, called a RVE, the size of which should be sufficiently 
large so that (1) it accommodates the unstable mode (2) the mechanical response is described well by the 
continuum constitutive relation. As long as in a given boundary value problem, we can identify such an 
RVE, the phonons can be used an indicator of material stability within the RVE. Once such an RVE is 
identified, the instability can be detected by phonon-based lattice instability analysis. 


Remark 2 : Since a soft mode phonon instability result from dynamical growth of certain modes of 
lattice-vibrations, a constrained boundary may suppress triggering of a short-wave instability in a nearby 
material element even if the local state of strain in the material element has reached the limiting condition. 
Therefore, the use of the limit surface in identification of an incipient failure is rigorously valid only in the 
interior, sufficiently far away from the boundary conditions. 


Remark 3 : An EE A simulation based on the continuum constitutive model of Sec. 0 collapses im¬ 
mediately after an elastic instability is encountered, since the hyperbolicity of the equations of motion no 
longer remains valid. However, since a continuum model is blind to short-wave instabilities, an EEA sim¬ 
ulation continues even if a short-wave instability is encountered, until it reaches the point of an elastic 
instability (see Eig. (13) in the following section). 


6.2. Prediction of incipient failure in blistering of a defect-free graphene sheet 

Employing the proposed scheme, we analyze the incipient failure of a defect-free graphene sheet in EEA 
simulations of idealized bulge-type tests. A bulge test involves blistering of a graphene sheet clamped at its 
boundary on top of a microcavity containing a pressurized gas. Two different geometries for the constrained 
boundary are considered in our work: elliptical and circular, while the characteristic dimension C, = 
remains same, where a and b are the sheet dimensions along the zigzag and armchair directions, respectively. 
Eor elliptical geometry, we choose a = 1.0 jam and 6 = 0.5 jam for the first configuration, and a = 0.5 jam and 
b = 1.0 jam for the second configuration; and for circular geometry, we choose a = b = l/\/2 /am. 



Eigure 11: Schematic of blistering (a) of an elliptic membrane and (b) of a circular membrane. The membrane is 
clamped at the boundary, and is subjected to a pressure (denoted by arrows) from below, resulting in the bulged 
configuration as shown in the figure. The state of strain at the center of the membrane is equi-biaxial in case of a 
circular membrane, whereas in case of elliptic membranes, it is biaxial, with larger strain in the minor direction. 
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Our objective is to determine the macroscopic and microscopic conditions at which the graphene sheet 
fails. Numerically implementing the constitutive response and the failure function in a commercial FEA 
program abaqus^ we explicitly determine the critical pressure-difference APq and the critical central de¬ 
flection the location of unstable material element; the local stress and strain at the unstable element; 
and the nature of instability. We demonstrate this procedure systematically for the blistering of the circular 
and elliptical graphene sheets. 

(a) Location of unstable material element and geometry of failed region 

To locate the unstable material element, we monitor the evolution of the limit function T = 7 ^(£:,e( 0 ))- 
7 over the suspended graphene sheet, as a function of externally-applied pressure-difference (shown in 
Fig. Q). The material element at which the limit function first vanishes denotes the unstable region, 
which will subsequently fail. This unstable region determines the geometry of the failed region. For exam¬ 
ple, in blistering of a circular graphene sheet — based on the negativity of the unstable region — we would 
expect that the fracture initiation occurs at the center in the form of a nearly circular void. 
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Figure 12 : Sequence of images showing the evolution of the limit function T = 7 ^(( 5 , 6 >)) — 7 with pressure- 
difference APq across the circular membrane. The central region with negative value of the limit function 
indicates the material elements that have reached the lattice-stability limit, which are at the verge of 
mechanical failure. 

(b) Pressure-difference and central deflection at the onset of instability 

The critical pressure-difference (APq ) and the critical central deflection {6^) are determined by locating 
the onset of instability on the pressure-deflection response curve (shown in Fig. ([T^). It should be noticed 
that the hyperbolicity of elastic-wave equations is essential for a finite element analysis step to continue; if 
an elastic instability is first reached along a loading path, as in the case of elliptical graphene sheet, then the 
condition of hyperbolicity ceases to hold, and the FEA step immediately terminates. On the other hand, if a 
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soft-mode instability is first reached on the loading path, the elastic-wave equation remains hyperbolic, and 
the FE steps continues until an elastic-instability along the loading path is reached. For example, during the 
blistering of a circular sheet, the microscopic soft-mode instability takes place prior to the elastic instability; 
the simulation continues past the soft-mode instability, terminating at the point of elastic-instability. We 
have indicated the initial encounter with both the soft-mode and the elastic instability on the pressure- 
deflection response curve. From Fig. (13), it clearly emerges that in the absence of the limit surface, an 
FEA simulation is bound to overestimate the mechanical strength of a crystal when the strength-limiting 
mechanism along the loading path is a short-wave instability. 



Figure 13: Pressure vs deflection response curves obtained from the elliptical and the circular blister simu¬ 
lations. On these curves, the long-wave (elastic) instabilities are located by a x; while the location of the 
short-wave (soft-mode) instability is denoted by +. 


(c) Local stress/strain at the unstable material element 

The local principal strains at the critical material element immediately after the onset of instability, 
and , are obtained by tracing the trajectory in the strain space, and obtaining the intersection of this 
trajectory with the limit surface, as shown in Fig. (14). Similarly, the local principal stresses at the critical 
element at the onset of instability, iS^in and are obtained by tracing the trajectory and locating its 

intersection with the limit surface in the stress space. These values are tabulated in Tab. (El). 


Geometry of constrained boundary 

‘5max(N/m) 

pc 

‘^max 

Nature of instability 

AP§ (MPa) 

Elliptical (a > b) 

37.74 

0.20 

Elastic 

72 

Elliptical (a < b) 

34.31 

0.175 

Elastic 

64 

Circular (a = b) 

31.28 

0.1475 

Soft-mode 

62 


Table 4: Local major principal stress and local major principal strain at the unstable material element as 
calculated by finite element analysis. It is noted that the strength is determined by a complex interplay 
between anisotropy, the geometry of the boundary condition and the nature of instability governing the 
failure. 

From Tab. 0, we note that the limiting stress and strain in the blistering simulations are inherently 
dependent on the geometry of the constrained boundary. Due to the fact that for certain geometries, e.g.. 
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when a ^ the failure is mediated by a microscopic soft-mode instability, the sensitivity of the failure 
stress/strain on the geometry of the constrained boundary can not be explained by the material anisotropy 
or the nonlinear softening in the material response alone. Thus, in a global sense, the proposed failure model 
is a useful method for assessment of correlation between the failure stress/strain and extrinsic factors. 

(d) Nature of instability: long-wave versus short-wave instabilities 

The proposed limit function captures lattice instabilities of all kinds; however, without resolving whether 
it is macroscopic (elastic) or microscopic (soft mode) in nature. By combining the proposed limit function 
with the acoustic tensor analysis, the distinction between macroscopic and microscopic instabilities can be 
readily made. To this end, we use the fact that an acoustic-tensor analysis captures macroscopic instabilities 
but not the microscopic ones [HI]. Thus, if the instability is indeed a macroscopic instability, the predicted 
onsets of instability from the failure function and the acoustic-tensor analysis coincide; however, when the 
instability is a soft mode, the predicted onset of failure from the limit-function precedes that from the 
acoustic-tensor analysis (see Fig.([l4|)). 



Figure 14: (a) Shown are the deformation paths followed by the critical material elements during the 

elliptical and the circular blister simulations, along with the lattice-stability limit surface. The intersection 
of the deformation path with the limit surface indicates the onset of failure, (b) Deformation paths for 
the elliptical (a < b) and the circular bulges, superposed on the ^ = 0 radial section of the limit surface. 
The long-wave (elastic) instabilities are located by a x; while the location of the short-wave (soft-mode) 
instability is denoted by +. For the elliptical membrane, the points of intersection of the trajectory with the 
limit surface and the elastic instability curve coincide, indicating that the failure is of elastic nature; while 
for the circular membrane, the intersection with the limit surface precedes that with the elastic instability 
curve, indicating that the failure is triggered by a soft mode instability. 
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7. Conclusion and discussion 


Onset of inelasticity or ‘failure’ in defect-free graphene crystals ensues from the deformation-induced in¬ 
stabilities of the underlying lattice. Though lattice-instabilities of all kinds can be assessed via an elaborate 
atomistic-level lattice-dynamical simulation, there remains a lack of a general continuum criterion capable 
of describing the conditions for the onset of inelasticity along an arbitrary deformation path, parametrized 
in terms of stress or strain. In addition to non-linearity and anisotropy in material response at large- 
deformations, a major difficulty in this regard arises due to the fact that different lattice-instability mecha¬ 
nisms govern the failure in along different loading paths. 

For a majority of deformation paths, material failure in crystals arises from a long-wave or elastic instability. 
Such elastic instabilities can easily be predicted by a continuum acoustic tensor analysis once an appropriate 
hyperelastic constitutive relation for the mechanical response of the material has been formulated [2I1I20]. 
Further, an elastic instability is readily seen in an FEA simulation, since upon reaching an elastic instabil¬ 
ity, the hyperbolicity of equations of motion ceases to hold and the unstable material element undergoes 
excessive distortion, annihilating the load carrying capacity of the material. However, there also remains 
a small but finite domain of deformation paths along which a short-wave (also called soft mode) emerges 
before an elastic instability is reached, and thus limiting the strength along those deformation paths m- 
Unfortunately, there exists no equivalent continuum criterion that could describe the onset of short-wave 
instabilities, and therefore a continuum FEA simulation remains blind to such instabilities. Thus, in the 
absence of a multi-scale resolution, an EEA simulation is bound to overestimate the mechanical strength of 
a crystal when the strength-limiting mechanism along the loading path is a short-wave instability. 

Here, we have presented for the first time the notion of a comprehensive lattice-stability limit surface — a 
continuum parametrization which describes the loss of internal stability of a defect-free graphene sheet — 
both in terms of stress and strain — with respect to perturbations of all kinds. Eurther, on this limit surface, 
we have identified the distinct regions that correspond to fundamentally different kinds of lattice-instability 
mechanisms. Such a continuum parametrization is valuable for several reasons: 

• Eirst, it duly accounts for all the complexities in crystal response, such as non-linearity, anisotropy and 
multiplicity of instability mechanisms, via one general continuum criterion: an instability appears when 
the magnitude of deviatoric strain reaches a critical value that depends upon the mean hydrostatic 
strain and the directionality of stretch; and the nature of this instability depends upon the mean 
hydrostatic strain level. 

• Secondly, despite being based on a continuum formulation, it possesses a multi-scale resolution. It can 
address the incipient failure of defect-free graphene crystals under an arbitrary loading condition and 
at realistic length scales — by the virtue of its continuum description — while the built-in microscopic 
lattice dynamical theory of stability analysis allows us to capture instabilities as-short-as the size of 
few unit-cells. 

• It is based on a novel interpolation scheme based on the symmetry-invariants of the logarithmic strain 
measure, which reduces the sampled deformation paths required for representation to a bare minimum 
of two: one along the zigzag direction and the other along the armchair direction. 

• Owing to its continuum representation, the limiting criterion is easily implemented in a continuum EEA 
simulation. Such an implementation is capable of performing an on-the-fly lattice stability analysis 
with an on-going finite element calculation, and can predict the location of the unstable material 
point, the load/displacement at the onset of instability, the local stress/strain, and the character of 
the instability, i.e., a short-wave or a long-wave in a realistic boundary-value problem. 

Employing the limit surface in an EEA simulation, we analyzed the failure of a defect-free suspended 
graphene sheet in the finite elements simulation of idealized bulge-type tests. Our simulations considered 
two different geometries of the constrained boundary, elliptical and circular, but with the same suspended 
area. We showed how different types of local lattice instabilities are reached in the respective boundary 
value problems under different levels of differential pressure and overall deflection. Erom these examples. 
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it emerged that, in the absence of a multi-scale resolution, a continuum-level finite element analysis is 
bound to overestimate the strength of graphene along certain deformation paths. Further, we showed that 
the proposed limit surface, owing to its continuum representation and multi-scale resolution, adequately 
remedies this deficiency. 

In a recent work [22], we have demonstrated the utility of the proposed limit surface in identifying a 
strain-shielding effect initiated by mechanically-activated covalent interactions at the graphene-indenter 
interface during nano-indentation of graphene. Employing the limit surface as the failure criterion in a FEA 
simulation, we explained how this strain-shielding effect can lead to a delay in the onset of instability in 
graphene during nano-indentation. 
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